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ABSTRACT 

We investigate the structure of self-gravitating polytropic stellar systems. We present a method which 
allows to obtain approximate analytical solutions, tp n + e (x) , of the nonlinear Poisson equation with the 
polytropic index n + e, given the solution ipnfe) with the polytropic index n, for any positive or negative 
e such that |e| <C 1. Application of this method to the spherically symmetric stellar polytropes with 
n ~ 5 yields the solutions which describe spatially bound systems if n < 5 and the formation of a second 
core if n > 5. A heuristic approximate expression for the radial profile is also presented. Due to the 
duality between stellar and gas polytropes, our results are valid for gaseous, self-gravitating polytropic 
systems (e.g., molecular clouds) with the index 7 ~ 6/5. Stability of such systems and observational 
consequences for both stellar and gaseous systems are discussed. 

Subject headings: galaxies: kinematics and dynamics — galaxies: star clusters — stars: kinematics — 
stars: formation — ISM: structure — ISM: clouds 



1. INTRODUCTION 

Galaxies, star clusters, and other stellar systems are sta- 
tionary assemblages of stars bound by a self-consistent, 
mean gravitational potential. Their structure is deter- 
mined by the nonlinear Poisson equation, which, in most 
cases, may be solved only numerically. Therefore, any an- 
alytical study acquires special interest. 

The simplest model of a gravitationally bound stellar 
system is a spherically symmetric polytrope with a power- 
law distribution function of stars. The gravitational po- 
tential of a polytropic sphere is determined by the Lane- 
Emden equation: 
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where r is the dimensionless spherical radius, ip(r) is the 
dimcnsionlcss potential, an d n is the polytropic index (see 



a finite radius — the star surface. Interestingly, the Lane- 
Emden equation appears also as an evolution equation for 
a scalar field in some ve rsions of the quintessence model 
( Liddlc fc Scherrer, 1999 ). Hence the present study is also 
of practical interest. 

There are two special cases for which nonsingular ana- 
lytical solutions of equation ([!]) are known, namely n = 1 
and n = 5. We performed a study of (|j) using a powerful 
mathematical techn ique — the Lie g roup analysis of dif- 
ferential equations ( Ibragimov, 19851) ■ This cumbersome 
analysis, the details of which are omitted from this publi- 
cation, yields that, for any n except n = 1 and 5, equation 
(0) is non-integrable in a closed form. 3 When n — 1, equa- 
tion (0) becomes the linear Helmholtz equ ation, and when 
n = 5, the solution has been discovered by Schuster (1883 ) 
and takes a simple form: 
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Tremaine, 1987 for more discussion). The den- 
sity is found from p = c n tp n , where c„ is a constant which 
depends on n. 

The stellar (collisionless) polytropes are related to hy- 
drostatic equilibria of self-gravitating spheres of a poly- 
tropic gas with the equation of state P = Kp' 1 , where 
7 = 1 + — , P is the gas pressure, and K is a, constant. The 
gas polytropes with 7 < 6/5 (i.e., n > 5) are usually used 
to model the structure of molecular clouds which has a ma- 
jor effect on the star formation that occurs within them. 
Note that the 7 — > 1 (n — > 00) limit corresponds to the 
case of an isothermal gas. The gravitationally stable poly- 
tropes 2 with 7 > 6/5 (i.e., n < 5) often serve as models of 
stars, because they can have a zero pressure boundary at 



1>(r) = [l + -r 2 



-1/2 



(2) 



For other values of n, the general solution may be found 
only numerically. 

The present study is motivated by a very peculiar depen- 
dence of the radial profiles of the gravitational potential 
and density on the value of n when n ~ 5. Several solu- 
tions of equation (|l|) with the "core" boundary conditions 
(i.e., -0(0) = 1, d r tp(0) = 0, where d r = £-) are repre- 
sented in Figure |l} A small change in n results in a drastic 
modification of the profiles. Solutions with n > 5 exhibit 
a multiple core structure, — the fact which has probably 
been known but never acknowledged in the literature. In 
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however, the absence of a Lie algebra often means nonintegrability of a differential equation. 
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this solutions, the regions where ip falls as a power-law 
are interchanged with flat cores where ip ~ constant. For 
n < 5, the solutions are spatially bound, i.e., ip and, hence, 
p vanish at a finite radius. The radius at which the second 
core begins to form in the n > 5 solutions is comparable 
to the radius at which the density vanishes in the n < 5 
solutions. 

In this paper we develop a functional series expansion 
technique which allows the derivation of solutions with the 
polytropic index n± |e|, given the solution with the index 
n (for any n and e, |e| <C 1), to any accuracy in powers of e 
and present it in §0. We apply this technique to the n = 5 
solution (||) and analytically derive a first-order (oc e 1 ) ap- 
proximate solution in §|3| This solution perfectly agrees 
with the exact numerical n < 5 solutions everywhere, and 
with the n > 5 solutions at radii smaller than the second 
core radius. We also provide a useful approximation for 
ip(r) for n ~ 5 to any power in e. These results are directly 
applicable to the polytropic gas systems with 7 ~ 6/5. In 
§[| we discuss stability of stellar polytropes. We discuss 
the results and make some observational predictions in 

2. ANALYTICAL STUDY OF POLYTROPIC SOLUTIONS: 
FORMALISM 

Let us rewrite equation (|l]) in a form of the nonlinear 
Poisson equation with some boundary conditions, 



V 2 i(x)+C(x) 
?Mx)| s = /(x), 



where V 2 is the three-dimensional Laplace operator, /(x) 
is some specified function, and S is a two-dimensional sur- 
face. The subscript n explicitly shows that ip n {~x) is a 
function of the parameter n, as well as three-dimensional 
coordinates x. Our goal is to find a solution tp n+e (x), given 
the solution VVi(x) for positive or negative e. For |e| <C 1, 
such a solution may be expressed as a functional expansion 
in powers of e: 
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Note, the derivative terms above are functions of coordi- 
nates x evaluated at a particular value of the polytropic 
index. Hereafter we use the short-hand notation for the 



derivatives, 



a">i(x) 



di" 



The 



= 0™^ n (x), m = 1,2,. 

derivatives d™ip n are found by solving equations (j^) and 
(||) below. Differentiating (||) with respect to n, we obtain 



V 2 (9 n ^ n )+nC _1 (9n^ 
d n tp n (-x.)\ s = 0, 



-C ln Vv, 



(5) 



where the last equation means that boundary conditions 
are independent of n. Equation (|^) is linear in the func- 
tion cLipn and, hence, is easier to solve than the origi- 
nal (||). Once the solution of a homogeneous equation 
[i.e., equation (|J) with vanishing right hand side] is found, 
it is straightforward to find a solution of inhomogeneous 
equation rtq), as shown in the Appendix. The quantity 
d„ip n represents the first-order correction in equation (|J). 
Higher-order terms are obtained by taking higher deriva- 
tives of (|^) and solving for appropriate functions. For in- 
stance, upon differentiating (|) m times, the equation for 



the m-th derivative, d™ip n , takes the form 



— — H n , 
Wn(x)U = 0. 



(ip n ,d n <ip n ,...,d™ Vn) 



(0) 



Here H n , m {- • •) is a function which depends on n, m, tp n , 
and all previously found functions d^ip n , 1 < k < m — 1. 
Note that the homogeneous solutions of equation (||) for 
all m are the same and identical to the homogeneous so- 
lution of equation (||). Hence, the calculation of d™ip n is 
staightforward (see Appendix). 

3. THE STRUCTURE OF STELLAR POLYTROPES WITH 

JV = 5±]e| 

Equations (||) , (||) allow the solutions ip n +e to be calcu- 
lated from a solution tp n . As mentioned above, there are 
two cases in which ip n is available. For n — 1, equation 
(||) is linear and the complete three-dimensional solution, 
V>i(x), is known. For n = 5, only the spherically symmet- 
ric solution, ipsir), has been discovered, which is given by 
equation (|2|) . Here we consider the latter case and find an 
analytical expression for the first-order correction. 

The function ip${r) is given by equation (|^). Equation 
(||) for n — 5 takes the form, 
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where we again introduced a short-hand notation (p{r) 



(7) 



d n ip5 



di 



i=5 



The general solution of this equation 



may be obtained if a special solution of the homogeneous 
equation [i.e., equation (pi) with the zero right hand side] is 
known, as shown in the Appendix. Since ip§ is a power-law 
of a rational function, we guess that a special solution of 
the homogeneous equation is a function of the same kind. 
We have the first homogeneous solution, 
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(8) 



The s eco nd homogeneous solution is obtained from equa- 
tion (Af) where A = 1/r 2 , as follows from equation (A3) 
It reads 
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(9) 



The general solution is found from equation (A2) by 
straightforward integration. The constants of integra- 
tion are found from boundary conditions. First, we re- 
quire the solution to be non-diverging everywhere; hence 
C'2 — 0. Second, we consider the solutions with a core, i.e., 
^ 5 (0) = V5+e(0) = 1; hence, 0(0) = d n ip 5 (0) = 0. The last 
condition yields C\ — \/3/12. Finally, we have 



-05+eM ~ ip 5 (r) + ed n tp 5 (r), 



where 



2\ -1/2 



(10) 



(11) 



JnMr)- 16r(r2 + 3 )3 /2 arc tan ^ 
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The approximate solution (|T^) is shown in Figure [l] by dot- 
ted curves for several values of e. The function d n yj5 is also 
shown by the dashed line. It is seen that solutions with 



n < 5 are very well approximated by the first-order expan- 



sion (ID) at all radii. In contrast, tor solutions with n > b, 
this approximation is valid for r less than or comparable 
to the second core radius, i.e., where the second "plateau" 
ends. To obtain the profiles at larger radii, higher-order 
terms in e must be kept. The next, e 2 -term can also be 
calculated analytically. It is, however, cumbersome and 
involves generalized hypergeometric functions and, there- 
fore, is not presented here. As for the higher-order cor- 
rections, the complexity of the resulting integrals seems to 
preclude further analytical treatment. 

We can now estimate the inner radius of the second 
core as follows. The asymptotic value of d n ip5 at r 3> 1 is 
d n ips(oo) = ir/32. The transition from a power-law decay 
to a flat core occurs when the first and second terms are 
of comparable magnitude, ■05( r *) — 7r|e|/32. Thus, the 
transition radius is 
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32^/3 



(13) 



Note that for e < this radius determines the size of the 
system, i.e., the radius at which the density vanishes. 

Without proof we suggest that the inner and outer core 
radii of other cores are also inversely proportional to the 
powers of e. Using this conjecture, we present a heuristic 
expression which approximates the n — 5 + |e| solution 
quite well. The expression is best written as a continuous 
fraction, 
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where r\ , ra , r% , . . . are constants independent of e which 
are to be found by fitting to the exact numerical solution 
with e — > 0. For the first two radii, we use the analytical 
expressions: r\ = y/3 and T2 = r*|e|. Truncating the frac- 
tion at the 77 term and performing the nonlinear fit to the 
exact solution for e = .001, we have 



V3 ~ 1.73, 
32V3/tt ~ 17.6, 
1.95 x 10 2 , 
9.29 x 10 2 , 
4.75 x 10 3 , 
1.60 x 10 4 , 
5.00 x 10 5 . 



(15) 



The polytropes with n = 5 — |e| are well described by the 
following approximate expression 



1 



y/l + r 2 /3 
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(16) 



which is easily obtained from equation (|10|). No higher 
order corrections are needed in this case. 

How do the n > 5 stellar polytropic spheres look to an 
observer? Assuming that the luminosity density is propor 
tional to the number density of stars , the surface bright- 
ness p rofile is calculated as follows (Binncy & Trcmainc 
"(198Tb ), 



I(R) = Jo / 
Jr 



p(r) r dr 
Vr 2 - R 2 ' 



(17) 



where R is the radial distance from the center of an object 
on the plane of the sky and Iq is a normalization. Note, if 
p(r) falls too slowly and the integral fails to converge, an 
upper limit must be replaced with an effective tidal radius 
to truncate the system. Using approximations (|I^)-(f5|) 
and that p n (r) oc ?/>™(r), we numerically determine the 
surface brightness profiles, as presented in Figure |2|. As 
one can see, photometric observations will reveal a bright 
core surrounded by a much larger and considerably fainter 
shell or "halo" throughout which the surface brightness is 
almost constant. Other halos which are located at larger 
radii are perhaps too faint to be observable. 

We can also estimate the profile of the velocity disper- 



sion of stars a 2 



which is analogous to the 



sound speed in a gas, a 2 (r) oc P/p oc i/j n (r). Therefore, 
the velocity dispersion decreases with radius and exhibits 
the same multi-core structure. In particular, the faint halo 
is cooler than the central core and looks nearly isothermal 
because a 2 ~ constant thoughout. 

4. STABILITY OF THE POLYTROPES 

The stability of self-gravitating polytropes depends on 
how the system responses to an adiabatic perturbation, 
i.e., whether the internal heat transfer through the system 
is long or short comp ared to the dynamical scale (McKee 
& Holliman II, 1999J). In the first case — a locally adia- 



batic polytrope — the polytropic index, which determines 
the global structure, may differ from the adiabatic index 
7ad, which determines thermodynamical properties. The 
second case represents a globally adiabatic polytrope in 
which 7 ad = 7. 

The stability analysis of nonsingular polytropic spheres 
usually requires numerical calculations, even for the sim- 
plest cases ( e.g., n = 5). Here is a brief s ummary of the 
results (see McKee fc Holliman II, 1999 and references 
therein). First, polytropes with 7 > 4/3 (i.e., n < 3) are 
unconditionally stable for any value of 7 a( j- Second, lo- 
cally adiabatic systems with 6/5 < 7 < 4/3 are stable if 
7ad > 4/3. A system with smaller 7 a( j, e.g., a globally adia- 
batic polytrope, is stable only if it is truncated at a certain 
radius and confined by an appropriate external pressure. 
Third, for the systems with 7 < 6/5, the truncation ra- 
dius depends on the value of 7 a d- This radius approaches 
infinity for 



7ad,i 



32 7 (2 - 7) 
(6- 7 ) 2 



< 7- 



(18) 



As one can see, globally adiabatic polytropes with 7 < 4/3 
are stable only if they are pressure confined (in an analogy 
with Bonnor-Ebcrt isothermal sphere). The truncation ra- 
dius follows from the solution of equation (54) of McKee 
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& Holliman II (1999) which in our dimensionless units 
reads as 

r r ^ / 4^ \ I _3 i 

m cr = / A-Kr 2 p{r)dr = f — ■— J ^ cr [p(r cr )] 2 " 2 . 

(19) 

Here factor determined numerically; it is close to 

unity for large n and monotonically approaches ~ 4.6 as 
n —>■ 3. For n = 5 and r cr ^> ri this equation yields 

Tor = — yff ■ (20) 

Mcr V 27T 

Comparing this to equation (|l3|), we see that r cr ~ r* 
corresponds to 

^fvf^^ 1 ' (21) 

i.e., n ~ 6. Therefore, the multiple core structure in a 
globally adiabatic system may be seen if its polytropic 
index is n ;> 6. The stability analysis of these systems 
requires the use of the density profiles given by equations 
( flCi| ) or ([l4]). The accuracy of our approximate solution 
degrades for e > 1. The qualitative behavior, however, 
remains correct (note that in the limit of an isothermal 
sphere e — > 00, some oscillations in the profile persist). On 
the other hand, the multiple core structure may be seen in 
pressure confined, locally adiabatic (7 a d 7^ 7) polytropic 
systems with n > 5 and 7 a d being sufficiently close to 
unity, so that r cr 3> r* . 

5. DISCUSSION 

In this paper we have demonstrated that the radial 
structure of polytropic self-gravitating spheres is very sen- 
sitive to the value of the polytropic index n when n ~ 5 
and have investigated it analytically. We developed a 
method which allows one to obtain approximate solutions, 
ipn+ei of the Lane-Emdcn equation, provided the solu- 
tion ip n is known and e is a positive or negative constant, 
|e| <C 1. We applied this method to n = 5 solution and ob- 
tained the first-order (in e) correction analytically. The ob- 
tained solution describes finite-size systems with n < 5 and 
the formation of the second core in systems with n > 5. 



Both effects occur at the radius r* ~ 32v3/7r|e|. We also 
provided a continuous fraction approximation for the poly- 
tropes with n > 5, equation (|l4|), and a simple approxi- 
mate formula for those with n < 5, equation (|l6|). 

We make the observational predictions by calculating 
numerically the surface brightness of a stellar polytropic 
system. A galaxy or a star cluster which is described by a 
polytrope with n ^> 5 will be observed as a bright core sur- 
rounded by a very extended faint halo or shell. This halo is 
cooler than the central core and looks nearly isothermal. 
Because of the large spatial size and low surface bright- 
ness this halo may easily be confused with a background. 
A way to disentangle the two is to probe a correlation be- 
tween the core radii, ri, r%, and (if detected) r% and the 
surface brightness distribution. 

Due to the duality between the stellar collisionless poly- 
tropes and the gas polytropes, our results are immedi- 
ately applicable to the latter with the adiabatic index 
7 = 1 + i ~ |. Such polytropic gas spheres are used 
in studies of molecular clouds and star formation. The so- 
lutions with 7 < 6/5 describe extended molecular clouds 
bounded by external pressure. Equation ([l4]) describes the 
multi-core structure of a cloud with 7 < 6/5. The den- 
sity and, hence, pressure is almost constant in core/halo 
regions. Therefore, if an external pressure, p ex t, is compa- 
rable to the pressure in one of the halos, then even a small 
variation in p ext will drastically change the radius of the 
cloud. Thus, a large dispersion of masses and sizes of the 
clouds embedded into nearly the same environments may 
indicate that the gas in the cloud has 7 <; 6/5, although 
other explanations are also possible. The gravitationally 
stable polytropes with 7 > 6/5 are often used as models 
of stars because the pressure vanishes at a finite radius, 
— the star's surface. If 7 p> 6/5, the radial profile is well 
described by equation ( |l6| ) and the radius of a star is given 
by equation (|l3|). 

We thank Ramesh Narayan for interesting discussions 
and an anonymous referee for very careful reading of the 
manuscript and many helpful suggestions. This work was 
supported in part by NASA grant NAG 5-2837 and NSF 
grant AST 9820686. 



APPENDIX 

GENERAL SOLUTION OF EQUATION (@) 
Let us consider the following differential equation, 

y" + f{x)y' + g{x)y = h{x), (Al) 

where prime denotes the x-derivative and fix), g(x), and h(x) are some functions. The general solution of this inhomo- 
geneous equation reads, 

y = Cm + C 2 yi + 2/2 J -gui dx-yi j ^y 2 dx, (A2) 

where C\ and C'2 are constants of integration, yi — y±(x) and y 2 = J/2 0*0 are t w0 linearly independent solutions of the 
homogeneous equation y" + f(x)y' + g{x)y = 0. Here also 

A = e-S f ^ dx = yi y 2 -y 2y 'i. (A3) 

4 There is a typo: in the first equality the power of c s is 3, not 3/2; cf., their equation (10). 
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If one nontrivial solution of the homogeneous equation, y\ , is known, then the second may be found from the equation 

A 



Wi = Vi 



dx. 



(A4) 
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Fig. 1. — Radial profiles of "05+e for several values of e. The solid lines are exact solutions of the Lane-Emden equation, the dotted lines 
are analytical approximate solutions, and the dashed line is the first-order term d n ip5(r). 

Fig. 2. — The surface brightness profile for n = 5.1, 5.01, and 5.001. 
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